Journal of Cheminformatics

A compound's toxicity can be gauged from the toxicities of the compounds it resembles. We showed that in 2003 for rat oral LD50: the compounds in a reference set of 13,645 that most resemble the query form a consortium and vote on its toxicity. This work returns to that method twenty-three years later, on public data, because the reference set used in 2003 is no longer public. Rebuilt on 11,989 curated rat oral LD50 values, the method answers for 67.0% of the held-out compounds carrying an exact value, and 92.1% of those answers land within one log unit; in 2003 the figures were 94% within one log unit on 57.3% of compounds. We then ran it over 28 further endpoints, each recording whether a compound showed a particular toxicity, from cytochrome inhibition to mutagenicity to clinical toxicity: 149,224 measurements on 34,162 compounds. The method answers only when a compound has close enough neighbors, and on those it correctly identifies a median 72.9% of the compounds that showed the toxicity and 96.5% of those that did not. Separately we asked whether a reverse screen carries information about these endpoints beyond molecular size and lipophilicity. We search a compound against a published index of 27,910 ligands solved inside proteins and count the protein families those sites belong to. A compound whose matches reach no protein family inhibits CYP3A4 2.5% of the time; one whose matches span sixteen or more families inhibits 51.2% of the time. Added to heavy atom count, molecular weight and cLogP, the screen raises the mean area under the ROC curve from 0.72 to 0.76 across the 28 endpoints. The OpenADMET consortium has since measured CYP3A4 inhibition on 2,335 compounds, and the same pattern holds. Compounds that resemble known binders of sixteen or more protein families inhibit CYP3A4 24.6% of the time; compounds whose matches reach no protein family inhibit it 7.3% of the time. The more protein families a compound resembles binders of, the more likely it is to inhibit CYP3A4. This association is strong enough to act on before anything is made. Scoring a library on how far its compounds reach across protein families, together with their size and lipophilicity, and setting aside the riskiest fifth removes two thirds of the androgen receptor binders and about half of the CYP3A4 inhibitors. Scoring needs only a two-dimensional structure, with no docking and no conformer generation at query time. Interactive scoring, the data and the pool are at toxpred.ai. Keywords: toxicity prediction; structural similarity; nearest neighbors; reverse screening; target breadth; promiscuity; LD50; cytochrome P450; virtual library triage

Deep learning models such as Chemprop have advanced quantitative molecular property prediction, but their reliance on large training sets limits use in data-scarce domains. We propose a framework that fine-tunes a general baseline model trained on publicly available data on small, class-specific datasets. The resulting models retain the baseline's generalization ability while gaining class-specific accuracy and produce probabilistic outputs that capture uncertainty in the training data. We demonstrate the approach on three toxicity classes defined by a common core structure, target, or mode of action: (i) organophosphates, (ii) androgen receptor antagonists, and (iii) estrogen receptor beta antagonists. Each fine-tuned model outperforms classical machine-learning methods and the EPA TEST tool. The probabilistic nature of the predictions enables prioritization of compounds for experimental validation and seamless integration with data streams of varying quality, supporting iterative decision-making in chemical safety and drug discovery.

BackgroundAdverse drug-drug interactions (DDIs) cause preventable hospitalizations, but exhaustive experimental screening of all drug pairs is infeasible. Many computational predictors rely on SMILES or other molecular representations, limiting their direct applicability to biologics and other non-small-molecule therapeutics. We present a structure-agnostic framework that combines semantic representations derived from Medical Subject Headings (MeSH) with graph-derived topology from a Drug-Target-Protein knowledge graph constructed from DrugBank and UniProt. We further investigate how variation in MeSH annotation depth affects predictive performance. ResultsDrugs are grouped according to their deepest MeSH annotation level (Low, Mid, or Deep), and performance is evaluated across the resulting interaction categories in transductive and inductive settings. The Intermediate ontology scope (Low+Mid) provides the most stable performance, while adding Deep-level terms offers limited and inconsistent benefit. Lightweight topological descriptors are integrated with MeSH features through instance-wise, dimension-specific latent-space gating, using curated reliable-negative pairs for supervision. Fusion improves mean performance over the MeSH-only baseline across all six categories in the transductive setting. Under induction, the clearest gains occur for Low-Low interactions ({Delta}AUROC = 0.056;{Delta} F1 = 0.137) and Low-Mid interactions ({Delta}AUROC = 0.077;{Delta} F1 = 0.114). ConclusionsMeSH annotation depth is associated with systematic variation in DDI prediction performance that aggregate evaluation can obscure. Graph-derived topology is particularly beneficial when ontology annotations are shallow. The framework provides a common, structure-agnostic representation compatible with both small-molecule and biologic therapeutics and supports first-pass DDI prioritization for subsequent expert assessment.

Machine learning (ML) models for molecular property prediction are increasingly deployed in drug discovery, yet their adoption in real-world scenarios requires an understanding of the conditions in which a model succeeds or fails. While standardized benchmarks are powerful instruments to measure and unlock progress in ML research, they should not be blindly treated as the end goal. Especially static and retrospective benchmarks, in which no true unknown test set is employed, limit our ability to robustly validate a model's performance. Building on the collective expertise of a cross-industry consortium, we present a model validation framework consisting of five recommendations that would enable the community to move beyond aggregate metrics toward understanding where and why molecular property prediction models fail. We connect evaluation choices to real-world applications and case studies encountered in pharmaceutical research. The framework proposes splitting strategies that mimic realistic distribution shifts and expose common failure modes. We apply the recommended framework to a recently released dataset of absorption, distribution, metabolism, and excretion (ADME) properties. Across two complementary model algorithms, our case studies reveal four distinct failure modes (extrapolation, interpolation, representation, and evaluation), showing that model errors arise not only from distribution shift but also from limitations in molecular representations. Our results show that commonly used evaluation protocols can significantly overestimate performance and may not detect important model failure modes. All software and data are released via https://github.com/srijitseal/polaris.

Accurately predicting cellular responses to drugs remains a challenge with the potential to reduce experimental screening costs and accelerate drug discovery. Current computational approaches represent compounds through chemical structures, which carry little information regarding their activity within biological systems. We show that gene expression responses, measured in a reference cell line, can instead serve as transferable representations of chemical perturbations. The proposed framework, BioPert, uses biological representations with a small neural network to predict transcriptional delta responses in cell lines of different lineages. Across the Tahoe-100M and LINCS L1000 datasets, BioPert outperforms molecular fingerprints and embeddings learned from chemical structure, which provide limited improvements over random controls. BioPert's prediction correlations reach 0.79 on Tahoe-100M, corresponding to an improvement of 0.34 over the next-best representation. On LINCS, performance varies with the experimental reproducibility of the test conditions, highlighting the impact of batch effects. Notably, BioPert considerably surpasses predictions obtained by copying the reference response, demonstrating that it captures context-dependent effects. Further analysis of a C32-cobimetinib case shows pathway-level accuracy of the predictions, including when the reference and target responses differ. These results open a new path for chemical perturbation prediction and could ultimately reduce the burden of phenotypic screening.

A three-dimensional pharmacophore fingerprint records the binding features a molecule can present. It is a description of a hand in search of a glove. Because it is defined by presented features instead of two-dimensional structure, it can identify pharmacophoric similarity between structurally distinct compounds, which is what scaffold hopping and non-obvious me-too design require. The descriptor has remained a niche tool because its cost is dominated by conformer generation. In the reference pipeline, generating 100 conformers requires 2.82 s of the 2.86 s needed to fingerprint one screening collection compound; the bit calculation requires 0.039 s. We therefore removed the conformational stage. PharmCast is a feedforward neural network that predicts all 10,549 bits of a PharmPrint ensemble fingerprint directly from a SMILES string. On the same machine, PharmCast generated pharmacophore fingerprints for two molecules and compared them in 0.584 ms, whereas the conventional conformer-based pipeline took 5.71 s. PharmCast version 10 was trained on 5,887,229 molecules drawn from a screening collection, activity-backed ChEMBL compounds from 142 to 1000 Da, and peptide loops excised from crystal structures. We evaluated 155,648 purchasable catalog compounds excluded from every training set, 139,700 activity-backed ChEMBL compounds not present in the version 10 training set, and 13,500 peptide loops reserved for testing. Median fingerprint error, Pearson r, and pairwise ranking accuracy were 0.008, 0.980, and 0.936 for screening collection chemistry; 0.016, 0.984, and 0.952 for loop peptides; and 0.027, 0.936, and 0.889 for activity-backed ChEMBL compounds. The reference calculation reproduces itself at an error of 0.006 and r of 0.995. Ensemble pharmacophore fingerprints can therefore be predicted from two-dimensional structure alone at a cost suitable for large-scale collection screening and virtual library exploration. Keywords: pharmacophore, fingerprint, scaffold hopping, surrogate model, virtual screening, applicability domain

Accurate ADMET prediction is essential for prioritizing compounds before costly experimental validation, yet ADMET tasks are highly heterogeneous. Properties such as solubility, permeability, protein binding, clearance, transporter activity and toxicity are governed by different molecular signals, ranging from local functional groups and physicochemical descriptors to bonded topology and three-dimensional geometry. Consequently, a single molecular representation or backbone is unlikely to be optimal across all ADMET tasks. We present Trimole-Hybrid, a task-wise multimodal framework that addresses ADMET heterogeneity by selecting or combining predictors built from complementary molecular representations. Trimole-Hybrid constructs a candidate pool of SMILES-, graph-, geometry-sensitive EPT/3D- and chemical descriptor-based predictors. For each task, Trimole-Hybrid selects the best-performing predictor to obtain the final prediction. On 22 Therapeutics Data Commons ADMET benchmarks, Trimole-Hybrid exceeded the public TDC top-1 methods on 10 tasks and ranked within the top 10 for 21 tasks. Ablation studies confirmed the contribution of both complementary multimodal molecular representations and task-specific ensemble strategies. In two small-molecule case studies, Trimole-Hybrid shows sensitivity to changes in essential functional motifs, suggesting its ability to capture ADMET-relevant molecular substructures.

High-quality molecular datasets are essential for reliable machine learning in cheminformatics and bioinformatics, yet dataset quality is rarely assessed systematically and its relationship with downstream model performance remains poorly understood. Here, we present MolJam, an open-source framework for quantitative assessment of molecular dataset quality across five dimensions-structural integrity, data quality, experimental information quality, chemical space coverage, and data distribution-using 12 standardized metrics. Application of MolJam to 11 MoleculeNet and eight ChEMBL-derived datasets revealed widespread and heterogeneous quality issues, including undefined stereochemistry in up to 70.72% of molecules, inconsistent molecular representations, and contradictory labels. We next asked whether improving these quality metrics necessarily improves machine learning performance. Refinement of the ESOL and Lipophilicity datasets increased their MolJam quality scores but produced mixed effects on predictive performance, suggesting a competing influence of reduced dataset size. Controlled ablation experiments further demonstrated that both dataset quality and data quantity contribute to model performance and, notably, that retaining molecules with incomplete stereochemical information can outperform their removal when the resulting gain in data quantity offsets the quality penalty. Thus, molecular dataset curation cannot be reduced to maximizing data cleanliness alone but requires balancing multiple dimensions of data quality against information loss. MolJam provides a standardized framework for diagnosing molecular dataset limitations, comparing benchmark quality, and quantitatively evaluating how data curation decisions influence downstream machine learning.

Protein-ligand interaction diagrams are a routine part of structural and medicinal chemistry, but the tools that produce them tend to force a choice: comprehensive detection with tabular output, publication-quality figures behind a licence, or a scripting environment that assumes expertise. PandaMap (Protein AND ligAnd interaction MAPper) is an open-source Python package that produces a 2D interaction diagram, an interactive 3D viewer, a text report, a machine-readable CSV, and a four-panel graphical summary from a single command. It reads PDB, mmCIF and PDBQT files, detects 15 interaction classes using crystallographically validated distance thresholds, and depends only on NumPy, Matplotlib, BioPython and Requests; RDKit improves the 2D ligand layout when present but is not required. Hydrogen bonds are filtered on the true D-H{middle dot} {middle dot} {middle dot} A angle when the structure contains explicit hydrogens, matching PLIPs 100{whitebullet} criterion on the same evidence, and on distance alone otherwise, with the provenance of each measurement recorded. We benchmarked the package on three complexes chosen for different chemistry: enolase with a phosphonate transition-state analogue (PDB 1ELS), the EGFR kinase with erlotinib (1M17), and aldose reductase with IDD594 (1US0). PandaMap recovers the contacts these structures are known for, including the EGFR hinge hydrogen bond to MET769 and the IDD594 bromine{middle dot} {middle dot} {middle dot} THR113 halogen bond, both at distances identical to PLIPs. All detection thresholds, scoring weights and the exact commands used are given in the Supplementary Information, and the release carries a regression suite covering each interaction class. PandaMap 4.3.0 is available on PyPI under the MIT licence.

Structure-based drug discovery is a widely used paradigm for the rational design of novel small molecule therapeutics. However, the benefits conferred by the use of structural information has seen limited adoption in machine learning, where ligand-only ("2D") models are still the industry standard for molecular property or binding affinity prediction. Structure-based ("3D") ML models for binding-affinity prediction promise to present a clear advantage, but have not yet overtaken existing 2D models. Here, we show that physics-based priors can improve predictive performance of structure-based models by comparing different model architectures with varying physical priors on several prediction tasks. We present the Modular Training and Evaluation of Neural Networks (mtenn) package, where we decompose affinity prediction into separate steps of embedding structure into learned representations and combining those embeddings into a predicted binding affinity. We consider both E(3)-invariant and E(3)-equivariant architectures to determine the importance of encoding roto-translational inductive biases, as well as different methods for combining learned embeddings. By first optimizing several aspects of model construction using the general purpose PDBBind dataset, we are able to improve the performance and data efficiency of structure-based models. When subsequently trained and evaluated on the COVID Moonshot small molecule drug discovery dataset, our tuned models perform on par with industry standard ligand-only models. Our decomposed model framework highlights that encoding some physical priors improves model performance, while more complex biases such as equivariance offer limited benefit. Additionally, structure-based models generalize better to an unseen target and display higher training efficiency. Overall, these results emphasize that structure-based models benefit from their ability to incorporate physics-informed constraints, giving promising directions for model architecture development. These results also suggest that the strength of these models may be in tasks specifically aimed at generalizability, providing guidelines for their use in early-stage drug discovery campaigns.

A compound in a cell meets a spectrum of proteins drawn from many families at once, while screening most often interrogates one target at a time. Two questions asked many times in a rank ordering workflow ultimately guide decisions on what gets made and what gets counter-screened, and both are comparative: which of two targets does a compound prefer, and which of two compounds does a target prefer. We built one model for each: the target comparison over a roster of 2,279 human proteins in 34 protein families, the compound comparison over 2,079 in 32. Each model is given two chemical structures and a sequence, or two sequences and a chemical structure, and returns which member of the pair is preferred together with how firmly it holds that view. No conformational analysis, protein structure, binding site or docked pose is used. Asked which of two targets a compound prefers, the model is correct 0.75 of the time over 8,689 held-out comparisons setting two families against each other, and 0.78 of the time over 32,738 comparisons between two targets of one family, rising to 0.94 and 0.95 on the most confident third of each. Asked which of two compounds a single target prefers, it is correct 0.71 of the time over 65,725 held-out comparisons, rising to 0.96 on the most confidently held. Neither compound in any of those comparisons appeared anywhere in training. Accuracy follows the gap between the two measurements. Where they sit within half a log unit the models are right 0.57 to 0.64 of the time, and where they differ by more than two logs, 0.90 to 0.92. The compound result holds across 31 protein families, not only the best-measured one. Within the chemistry and the targets they were built on, these models rank compounds and rank targets well. Both models can be explored and downloaded at familyfoundationmodel.com. Keywords: target preference; compound preference; pairwise comparison; polypharmacology; off-target triage; ESM2; random forest; structure-free prediction; ChEMBL

Protein-cofactor selectivity prediction is important for enzyme annotation, metabolic engineering, and biocatalyst discovery, but it remains challenging because experimentally supported associations are limited and unevenly distributed, and because biologically relevant differences can be subtler than those captured by general protein-ligand interaction models. Here, we present MINT (Multimodal Interaction Network with Teacher supervision), a multimodal framework for cofactor specificity prediction. MINT combines protein sequence features with cofactor sequence and graph representations and is trained with a gradual adaptation strategy that transfers useful information from a pretrained protein-ligand interaction model to the cofactor prediction task. On a fixed balanced held-out test set, MINT achieved the best overall performance among all compared methods (AUC = 0.932), and this ranking was reproduced in an independent family-stratified evaluation at the natural positive rate (pooled AUC = 0.939). Further analyses suggest that the pretrained model provides transferable structural information, while the gradual training scheme helps adapt that information to cofactor specificity. In a case study of the chemically near-identical NAD/NADP pair, MINT recovered 95.0-98.6% of true binders and resolved a single-phosphate-level distinction with direct relevance to engineering enzyme cofactor specificity. These results show that gradual transfer learning is an effective approach to cofactor specificity prediction when annotation is limited.

Molecular pretraining offers a route to learning transferable chemical representations from unlabelled data. However, existing approaches pretrained on isolated molecular structures struggle to generalize to reaction-specific tasks because their pretraining objectives provide limited information about chemical transformations. A key challenge is to integrate this reaction information into chemical representations. Here we introduce LARK, a unified framework for learning atom-molecule-reaction knowledge from chemical transformations. Using approximately 0.6 million reactions encompassing 0.88 million unique molecular components, LARK combines bidirectional bond-electron reconstruction with structural and geometric supervision through a role-aware hierarchical hypergraph Transformer that links local atom representations to reaction context. The learned representations transfer to both reaction prediction and single-molecule property prediction, with strong performance across 42 tasks covering reaction outcomes, physicochemical properties, toxicity, metabolism and pharmacokinetics. Reactioncentre analyses show that LARK identifies where bonds break and form, while targeted interventions reveal that the learned atom representations contribute to reconstructing these changes. These findings indicate that LARK learns chemically informative representations of substructures involved in molecular transformations. Together, the results establish chemical transformations as a source of structured supervision for representation learning across molecular and reaction scales.

Drug discovery (DD) is a complex, time-consuming, and resource-intensive process that involves the identification of therapeutic targets, selection of bioactive compounds, and extensive experimental validation. The discovery of promising therapeutic compounds from large libraries of phytochemicals and synthetic molecules remains a major challenge in modern drug development. Screening millions of compounds through conventional experimental approaches requires substantial time, cost, and computational resources. In recent years, in- silico molecular docking has emerged as an important computational approach for predicting interactions between small molecules and target proteins, thereby helping researchers prioritize promising compounds for further investigation. Several molecular docking webservers, including iScreen, SwissDock, CB-Dock2, DockThor, and MTiOpenScreen, have been developed to support virtual screening studies. However, many currently available platforms still face some important limitations. Most existing tools lack integrated repositories of medicinal plant-derived phytochemicals and organism-derived bioactive compounds, automated mapping between plants and their associated phytochemicals, and flexible ligand retrieval using chemical names, SMILES strings, PubChem CIDs, or drug names. In addition, many platforms require extensive manual protein and ligand preparation, provide limited support for AlphaFold-predicted protein structures, and lack efficient large-scale multi-target virtual screening. Most existing docking platforms offer limited support for interactive inspection of docked protein-ligand complexes, often requiring users to download the results and analyse them using external molecular visualization software. To address these limitations, we developed FlexAutoDock, an automated cloud-based molecular docking platform that provides a unified environment for protein-ligand docking and large-scale virtual screening. Unlike existing web servers, FlexAutoDock integrates curated repositories of medicinal plant- derived phytochemicals, organism-derived bioactive compounds, and synthetic compounds from the ZINC database while supporting flexible ligand acquisition through medicinal plant or organism selection, chemical names, SMILES strings, PubChem CIDs, and drug-name queries. The platform further streamlines the docking workflow through automated protein structure retrieval from the Protein Data Bank and AlphaFold databases, receptor and ligand preparation, chain-specific protein selection, blind and site-specific docking, interactive visualization of predicted protein-ligand complexes, and scalable multi-target virtual screening. The resulting platform enables rapid, flexible, and large-scale virtual screening while simplifying the molecular docking workflow, providing researchers with an accessible computational resource for accelerating early-stage drug discovery. FlexAutoDock offers a fast, reliable, and accessible computational platform for molecular docking and virtual screening, freely available to the scientific community at http://103.99.177.82:3000/.

Predicting ligand activity across G protein-coupled receptors (GPCRs) requires models that capture both molecular structure and receptor-specific information while remaining robust to chemical and target-domain shifts. We developed a multimodal graph neural network that combines explicit ligand molecular graphs with frozen protein language model representations of GPCR sequences and jointly predicts quantitative pActivity and binary activity. The model was trained on 271,739 curated ligand-GPCR pairs spanning 183,694 ligands and 216 human GPCRs and evaluated using random, Bemis-Murcko scaffold, and strict cold-ligand partitions. With ESM-2 650M receptor embeddings, the selected model achieved mean absolute errors of 0.513 {+/-} 0.006, 0.540 {+/-} 0.006, and 0.641 {+/-} 0.005 pActivity units under random, cold-ligand, and scaffold evaluation, respectively, substantially outperforming a protein-aware fixed-feature multilayer perceptron and consistently improving quantitative prediction over a matched GINE reference. Removing receptor embeddings markedly degraded both regression and classification, whereas differences among ESM-2 35M, ESM-2 650M, and ProtT5 were comparatively small. Independent evaluation on 6,319 ChEMBL 37/BindingDB pairs revealed a substantial external domain shift, with mean absolute error increasing to approximately 0.94-0.95 despite chemically stringent internal validation. Nevertheless, receptor-specific information remained highly informative in a DRD2-DRD3 selectivity analysis: the cold-ligand ensemble achieved R^2 = 0.853, Spearman {rho} = 0.916, and ROC-AUC = 0.966 for strong DRD3 selectivity, whereas the receptor-independent ablation approached chance-level discrimination. These results establish protein-language-model-conditioned molecular graph learning as a scalable strategy for GPCR bioactivity prediction while identifying chemotype novelty, external-domain transfer, and richer ligand-receptor interaction representations as key remaining challenges.

Accurate identification of drug-target interactions is consequential for novel drug discovery and development. Deep learning methods for drug-target affinity (DTA) prediction have shown great promise in accelerating drug discovery and reducing development costs. Although graph neural networks have improved drug representation learning for DTA prediction tasks, many models still struggle to effectively and efficiently capture protein information, limiting overall prediction accuracy. In this work, we systematically evaluate the impact of pre-trained protein language models (PLMs) on the downstream task of predicting binding affinity between drugs and target proteins. We design multiple experiments across four different molecular representation backbones and assess the effect of incorporating PLM embeddings, comparing their performance to classical 1D convolution methods. We evaluate four families of PLMs which we integrate into PLM-GraphDTA, each built on distinct architectures and optimized for different tasks, including structure prediction, function prediction, and sequence unmasking. Additionally, we evaluate DeepGraphDTA, an architectural modification of the baseline convolution method designed to improve protein representation learning. The models are evaluated on two benchmark datasets, Davis and KIBA, using concordance index (CI) and mean squared error (MSE) as performance metrics. We further evaluate the generalization power of each model using cold-start train and test splits, and analyze the per-protein contribution to total CI. The results indicate simple architectural modifications to traditional convolution methods may be sufficient to bridge the gap to large pre-trained PLMs.

Structure-based scoring functions leveraging machine learning have recently demonstrated superior performance over classical scoring functions, particularly on virtual screening benchmarks. However, due to the fundamental differences between their underlying model principles and architectures, it remains unclear to what extent performance stems from an understanding of molecular binding or from exploitation of systemic biases. Thus, disentangling the factors underlying benchmark performance is essential for determining whether a scoring function will generalize to novel chemical space and succeed in prospective drug discovery. To address this need, we present a case study investigating the nature and impact of systemic biases on benchmark comparisons between different scoring function paradigms. By systematically analyzing the evaluation workflows of prominent models, we reveal pocket bias, a form of spatial coordinate frame leakage arising from static binding pocket extraction, which artificially inflates benchmark performance. To progressively eliminate these sources of bias, we benchmarked two selected graph neural network scoring functions against two minmalist machine learning models and a classical scoring function under four increasingly stringent evaluation levels, successively removing pocket bias, reducing structural data leakage, and finally evaluating on out-of-distribution (OOD) protein targets. Upon removal of pocket bias and structural data leakage, the performance of all machine learning models dropped substantially. When evaluated on out-of-distribution protein families, the classical baseline AutoDock Vina outperformed the machine learning models in five of seven virtual screening tasks and dominated the docking power evaluation. Our findings indicate that benchmark performance can be heavily shaped by evaluation design and dataset artifacts, potentially overshadowing algorithmic improvements. While the tested machine learning models remain heavily dependent on encountering familiar data distributions to achieve competitive results, AutoDock Vina demonstrated superior generalization capacity on OOD targets. This work underscores the critical need for rigorous, artifact-free benchmarking protocols to guide the development of truly prospective machine learning models for virtual screening.

Identifying protein targets from phenotype-first or mechanism-uncertain compounds remains difficult because proteome-scale docking can generate thousands of structurally plausible interactions per compound. We developed a workflow that converts ranked proteome-scale docking profiles into stable Gene Ontology (GO) Biological Process enrichment signatures and evaluates whether those signatures can reduce the candidate target space while preferentially retaining known drug-target relationships. Docking targets were retained at the PDB-chain level, mapped to unique human gene identities, and analyzed with PANTHER overrepresentation against the screened structural gene universe. GO enrichment was evaluated from the top 25 through 505 ranked proteins in increments of 10, and a stable compound-level GO profile was selected using a Jaccard stability threshold of 0.80 across three consecutive transitions. Benchmarking used Yamanishi drug-target interactions, with 682 mapped compounds assigned to a prespecified development/held-out split (552/130) and evaluated at canonical, mechanistic, and fine-mechanistic biological resolutions. In the full corrected benchmark, 642 compounds with usable canonical GO profiles showed greater within-class than between-class similarity (0.1019 versus 0.0931; delta = 0.0088; 100,000-permutation p = 0.00222), and canonical class explained 1.34% of multivariate GO-profile variation by PERMANOVA (p < 1e-4). Fine-mechanistic labels showed stronger organization in the full dataset (delta = 0.0245; PERMANOVA R-squared = 0.0976; both p < 1e-4). Held-out validation was more modest and metric-dependent: mechanistic labels were significant by PERMANOVA (R-squared = 0.0626, p = 0.0437), whereas the frozen fine-mechanistic analysis showed greater within-class similarity (0.1350 versus 0.1110; p = 0.038) and significant nearest-neighbor recovery (p = 0.0495), but not significant PERMANOVA (p = 0.119). The principal held-out search-space experiment evaluated 107 compounds, 753,492 candidate protein rows, and 424 represented gold-standard targets. A direct GO gate retained 1.26% of candidates while retaining 11.32% of known targets (8.96-fold enrichment); ontology-propagated GO associations retained 5.10% of candidates and 21.46% of known targets (4.21-fold enrichment). At matched candidate-space sizes, GO-associated prioritization retained 27.59% versus 22.41% of known targets at approximately 5% of candidates, 39.39% versus 36.08% at 10%, and 58.73% versus 56.13% at 20%. By contrast, additive protein-level GO reranking was heterogeneous: among 605 evaluable compounds, 21.7% improved their best known-target rank, but mean reciprocal rank decreased from 0.0276 to 0.0189. These results support GO enrichment as an intermediate biological search-space reduction and prioritization layer rather than a universal direct target-scoring function. Keywords: proteome-scale docking; Gene Ontology; target discovery; targetome; PANTHER; target fishing; biological filtering; search-space reduction; Yamanishi benchmark

Large language model (LLM) agents are increasingly capable of carrying out autonomous computational research, but it remains unclear whether they can develop molecular modeling methods that compete with strong human-developed approaches. Here, we evaluate autonomous method development across four settings: Therapeutics Data Commons (TDC) ADMET tasks, the OpenADMET ExpansionRx Challenge, the activity prediction track of the OpenADMET PXR Induction Challenge, and the Practical Molecular Optimization (PMO) benchmark. Agents followed the prescribed data splits, metrics, and evaluation protocols for each benchmark. On TDC, the best agent-developed score across configurations improved on the leaderboard reference on 18 of 21 tasks and on the combined leaderboard and peer-reviewed reference on 17, while individual Codex configurations improved on the combined reference on 13 tasks. On ExpansionRx, five-agent and single-agent Codex configurations achieved overall scores corresponding to fourth and sixth place relative to the published final leaderboard, with the five-agent configuration reaching top-five performance on seven of nine endpoints. On PXR, five-agent Codex achieved a mean absolute error slightly outperforming the best published challenge entry. On PMO, five-agent Codex developed a molecular optimizer that outperformed existing methods evaluated under the benchmark's standard setup; other methods obtained higher reported scores under modified conditions, including additional molecular pretraining data or a different oracle-call budget. Together, these results show that current LLM agents can autonomously develop high-performing molecular prediction and optimization methods across substantially different drug-discovery settings, in several cases matching or exceeding leading human-developed approaches.

Within the ADMET continuum, metabolism and excretion (ME) form a critical bridge between drug exposure established by absorption and distribution and downstream efficacy and toxicity. However, no existing framework for drug ME prediction has simultaneously achieved broad coverage of endpoints and robust data recency and completeness. Here, we developed MetExPred, a multi-view prediction framework that comprises 17 classification endpoints and two regression endpoints, covering the overall drug ME process. The framework combines sequence-based and graph-based molecular representations, while optionally incorporating ESM-2 protein representations when experimentally annotated targets are available. A protein-aware masking strategy enables the same architecture to operate in both molecular-only and target-enhanced settings, and multi-view attention adaptively integrates the available representations. Across the classification tasks, MetExPred achieved mean AUROC, AUPRC and F1 scores of 0.844, 0.734 and 0.697, respectively. For clearance and half-life prediction, the model achieved RMSE values of 0.686 and 0.668. MetExPred showed the strongest average performance among the evaluated baselines, while ablation studies confirmed the complementary contributions of molecular sequence, graph and protein information. These results provide a unified and flexible modeling framework for systematic ME prediction of drug compounds, enabling early-stage virtual screening and pharmacokinetic assessment during lead optimization.